Section 9 Model Parameterization
These sections are not detailed enough to warrant their own dedicated chapter.
9.1 Parameter Changes
These sections mostly involve changing a few parameter values and are therefore in small sub sections without R-code. Some of this might change later.
9.1.1 Channel parameters revised
This an ongoing unresolved issue #49.
So far nothing has been changed.
9.1.2 Crop parameters verified
The crop database has been updated and stored in the following file:
"model_data/cs10_setup/swat_input_mod/plants.plt"Discussions to this topic can be found in issue #27.
From the whole database, we only use crops defined in our management files:
Crop management: oats, barl, wwht, swht, pota, fesc
Generic: frst, fesc and others which are minor
The parameters for these crops have been updated to reflect the colder growing conditions.
9.1.3 Soil physical parameters in final form
This has been mentioned in issue #28, and has been closed without discussion. Seems like the parameters are in final form, but no further documentation exists at this point.
9.1.4 Soil chemical parameters in final form
This step will be verified in issue #46. It has been discussed in #19.
The following changes have been implemented, but are subject to change:
| Parameter | Default | CS10 |
|---|---|---|
lab_p |
5 | 20 |
nitrate |
7 | 8.5 |
inorgp |
3.5 | 35.1 |
watersol_p |
0.15 | 0.4 |
nutrients_sol <- readLines("model_data/cs10_setup/swat_input/nutrients.sol")
header <- nutrients_sol[1]
nutrients_sol_df <- read.table("model_data/cs10_setup/swat_input/nutrients.sol", header = T, skip = 1, sep = "", fill = T, as.is = T, colClasses = "character")
nutrients_sol_df$lab_p <- 20
nutrients_sol_df$nitrate <- 8.5
nutrients_sol_df$inorgp <- 35.1
nutrients_sol_df$watersol_p <- 0.4
write.table(
nutrients_sol_df,
"model_data/cs10_setup/swat_input_mod/nutrients.sol",
quote = F,
sep = " ",
row.names = F,
)
header_new <- paste(header, "and modified by the CS10 workflow on", Sys.time())
nutrients_sol <- readLines("model_data/cs10_setup/swat_input_mod/nutrients.sol")
writeLines(c(header_new, nutrients_sol), "model_data/cs10_setup/swat_input_mod/nutrients.sol")Testing if the setup still works:
simout <- test_swat_mod()
cat(tail(simout), sep = "\n")9.1.7 Point sources parameters added
Is an ongoing issue in #21
9.1.8 Tile drainage parameters defined
Is an ongoing issue in #53.
old_path <- "model_data/cs10_setup/swat_input/tiledrain.str"
new_path <- "model_data/cs10_setup/swat_input_mod/tiledrain.str"
tiledrain_str <- readLines(old_path)
header <- tiledrain_str[1]
header <- paste(header, "then modified by CS10 workflow at", Sys.time())
tiledrain_df <-
read.table(
old_path,
sep = "",
header = T,
skip = 1,
fill = T,
colClasses = "character"
)
tiledrain_df$dp <- 800 # from 1000
tiledrain_df$t_fc <- 12 # from 24
tiledrain_df$lag <- 30 # from 96
tiledrain_df$rad <- 200 # from 100
tiledrain_df$dist <- 8000 # from 30
tiledrain_df$drain <- 40 # from 10
tiledrain_df$pump <- 0 # from 1
tiledrain_df$lat_ksat <- 2 # from 2 (no change)
write.table(header, file = new_path, quote = F, col.names = F, row.names = F)
write.table(
x = tiledrain_df,
file = new_path,
sep = " ",
quote = F,
append = T,
row.names = F
)## Warning in write.table(x = tiledrain_df, file = new_path, sep = " ", quote = F,
## : appending column names to fileTesting if the setup still works:
simout <- test_swat_mod()
cat(tail(simout), sep = "\n")9.1.10 Additional settings verified
Refers to chapter 3.11 in The Protocol and files parameters.bsn codes.bsn
Has been discussed in
#22
The following changes have been made:
| Parameter | Default | CS10 |
|---|---|---|
pet |
1 | 2 |
rte_cha |
0 | 1 |
cn |
0 | 1 |
tiledrain |
0 | 1 |
soil_p |
0 | 1 |
atmo_dep |
a | m |
old_path <- "model_data/cs10_setup/swat_input/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"
codes_bsn <- readLines(old_path)
header <- codes_bsn[1]
header <- paste(header, "then modified by CS10 workflow on", Sys.time())
codes_df <- read.table(old_path, skip = 1, header = T, sep = "",
colClasses = "character")PET method (pet) The OPTAIN project recommends the PET calculation of Pennman-Monteith, which is code 2.
codes_df$pet <- 2 Channel routing network (rte_cha) Recommended is to start with Muskingum (code 1) and apply variable storage method if it causes problems.
codes_df$rte_cha <- 1Stream water quality (wq_cha) Recommended to test both for OPTAIN (1/0)
Daily curve number calculation (cn) Code 1 is recommended. (Plant ET)
codes_df$cn <- 1OPTAIN recommends new version
codes_df$soil_p <- 1Lapse rate is being tested in REF
Plant growth stress is being testing in REF
Tile drainage
This should be set to 1, but that causes a crash. We will fix this bug in a dedicated secton (REF)
codes_df$tiledrain <- 0Atmosdep was done in REF and is annual
codes_df$atmo_dep <- "y"Writing the changes
write.table(header, file = new_path, col.names = F, row.names = F, quote = F)
write.table(
codes_df,
file = new_path,
col.names = T,
append = T,
quote = F,
sep = " ",
row.names = F
)## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to fileTesting if the setup still works:
simout <- test_swat_mod()
cat(tail(simout), sep = "\n")9.2 Fixing Tiledrain
run_this_chapter = FALSEProblem: tiledrain=1 in codes.bsn causes crash.
Reason: Some HRUs have drains, but are on a shallow soil type. This causes the drains to be below the soil, leading to a (nondescript) error.
Likely source of error: Misalignment of the soil and land use map.
Discussion of this issue can be found in #53
require(mapview)
require(sf)
require(stringr)
require(dplyr)
require(purrr)
require(DT)9.2.1 Finding the problem HRUs
Our drains are set to 80cm depth. Which of our soils are less than that?
usersoil <- readr::read_csv("model_data/input/soil/UserSoil_Krakstad.csv", show_col_types = F)
shallow_soils <- usersoil$SNAM[which(usersoil$SOL_ZMX<800)]
print(shallow_soils)## [1] "Mountain" "CMlen-dy" "Antropogenic" "Sea_dep_thin" "Humic"
## [6] "LVlen-sl" "End_moraine" "RGlen-dy-ar" "RGlep-dy" "STlen-dy"
## [11] "STlen-lv-sl" "STlen-rt-sl" "STlen-um-sl"Which of our HRUs use these soils, and are drained?
hru_data <- read.table( "model_data/cs10_setup/swat_input/hru-data.hru",
skip = 1, header = T, sep = "")
shallow_hrus <- which(hru_data$soil %in% shallow_soils)
drained_hrus <- grepl(x = hru_data$lu_mgt, pattern = "_drn") %>% which()
# Intersection of HRUs which are both shallow and drained
problem_hrus <- base::intersect(shallow_hrus, drained_hrus)
length(problem_hrus)## [1] 16Just 16 HRUs have an issue. Lets get some info on these 16:
problem_hru_data <- hru_data[problem_hrus,]
datatable(problem_hru_data)We need to link the soil type to the hru vector file
hru_shp <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
hru_soil <- hru_data %>% select(soil, name)
hru_shp <- left_join(hru_shp, hru_soil, by = "name")And separate out our problematic HRUs and affected LUs
problem_shps <- hru_shp %>% filter(name %in% problem_hru_data$name)
problem_hru_ids <- problem_hru_data$lu_mgt %>% str_remove("_drn_lum")
problem_lus <- hru_shp %>% filter(type %in% problem_hru_ids)Next we remove the problem HRUs and affect LUs from our hru_shp (this is for the upocoming map, to remove overlaps).
hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_lus$name) == FALSE)
hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_shps$name) == FALSE)We also remove the problem HRUs from the affected LUs
problem_lus <- problem_lus %>% filter((problem_lus$name %in% problem_shps$name) == FALSE)Data processing done, lets have a look:
hru_map <-mapview(
hru_shp,
map.types = "Esri.WorldImagery",
color = "green",
lwd = 1,
legend = FALSE,
zcol = "soil",
alpha.region = 0
)
problem_map <-
mapview(
problem_shps,
color = "red",
lwd = 3,
legend = FALSE,
zcol = "name",
layer.name = "Problem HRUs",
alpha.region = 0
)
problem_lus_map <-
mapview(
problem_lus,
color = "orange",
lwd = 1,
legend = FALSE,
zcol = "name",
layer.name = "Problem LUs",
alpha.region = 0
)
full_map <- hru_map + problem_lus_map+ problem_map
full_map